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The widely used density matrix renormalization group (DRMG) method often fails to converge in 
systems with multiple length scales, such as lattice discretizations of continuum models and dilute or 
weakly doped lattice models. The local optimization employed by DMRG to optimize the wave func- 
tion is ineffective in updating large-scale features. Here we present a multigrid algorithm that solves 
these convergence problems by optimizing the wave function at different spatial resolutions. We 
demonstrate its effectiveness by simulating bosons in continuous space, and study non-adiabaticity 
when ramping up the amplitude of an optical lattice. The algorithm can be generalized to tensor 
network methods, and be combined with the contractor renormalization group (CORE) method to 
study dilute and weakly doped lattice models. 

PACS numbers: 02.70.-c, 37.10.Jk, 67.85. Hj, 71.27.4-a 
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The optimization of variational wave functions is gen- 
erally a very difficult problem. In the specific case 
of matrix-product states (MPS) [:] the density-matrix 
renormalization group algorithm [2-5] often reliably and 
efficiently optimizes these wave functions to find a good 
approximation of the ground state. While most efficient 
in one dimension, it can be applied to medium-sized two- 
dimensional systems [ij] , and has been generalized to cal- 
culate time-dependent [~-'^^] and finite temperature prop- 
erties [10-12]. 

In systems with multiple length scales, however, the 
DMRG algorithm often fails to converge, as the local op- 
timizations that are at the core of DMRG are ineffective 
in optimizing large-scale features of the wave function. 
Especially in dilute systems where the inter-particle dis- 
tance is large compared to the lattice spacing the con- 
vergence of the density profile can be very slow. Sys- 
tems with multiple length scales suffering from this prob- 
lem arise from lattice discretizations of continuum mod- 
els [13], or in weakly doped lattice models where the hole 
density exhibits the same convergence problems. The 
first situation was recently discussed in Rcf. [ ]. 

Similar convergence problems are also known in other 
fields, e.g. when solving partial differential equations [io], 
lattice field theories [IG] or electronic structures [IT], 
and have there been overcome by multigrid approaches. 
Multigrid methods use a hierarchy of discretizations, as 
sketched in Fig. 1. Starting from the target problem on 
the finest grid (or a lattice model), the system is mapped 
to hierarchy of coarser grids. An approximate solution 
of the smallest problem on the coarsest grid is then used 
to initialize optimizations of the problem on the next 
finer grid and this process is iterated down to the finest 
grid. This method can substantially speed up a calcu- 
lation since the large scale features converge quickly on 
the coarsest grid, and the following calculations on finer 
grids only need to optimize local features at the scale of 
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Figure 1. Multigrid DMRG illustrated for bosons in an opti- 
cal lattice V{x). The DMRG algorithm converges fast for the 
rather dense system at the coarsest grid (shown on the top), 
and this solution is then used to iteratively initialize DMRG 
calculations on finer grids, which substantially speeds up con- 
vergence. The filling of the circles illustrates the probability 
for a particle to be at that site. 



the respective grid spacing. 

In this Letter we develop a multigrid DMRG (MG- 
DMRG) algorithm to solve the above-mentioned con- 
vergence problems in DMRG calculations. As a first 
application and demonstration of the effectiveness of 
the algorithm we present results for bosons in continu- 
ous space where MG-DMRG enables the study of non- 
adiabaticities when slowly ramping up the amplitude of 
an optical lattice. 

We start the description of the MG-DMRG algorithm 
by reviewing MPS wave functions on a chain of L sites; 
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used in DMRG. They are characterized by a polynomial 
number c>c LM^ of variational parameters, the M x M 
matrices A^^ . In one dimension a good approximation for 
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Figure 2. Density profiles of a continuum system of bosons 
in an optical lattice consisting of L = 32 unit cells. The top 
and middle panels shows the non-converged results obtained 
for A'^ = 32 grid points per unit cell after 12 sweeps of the 
DMRG algorithm with two different initial states: initial state 
1 is a random state, initial state 2 is obtained from an infinite- 
size growing procedure as implemented in the ALPS DMRG 
code [-']■ The bottom panel shows the multigrid result. 



low-energy states can be obtained by MPS wave functions 
with a fixed or at most polynomially growing M [Is, 19]. 

While optimizing MPS wave functions to obtain a vari- 
ational estimate for the ground state is a hard non-linear 
problem, the DMRG algorithm is very effective in many 
cases. It iteratively optimizes one or two of the matrices 
A"^' while keeping all other matrices fixed, and sweeps 
back and forth along a (quasi) one-dimensional system 
until convergence is achieved. For a recent review and im- 
plementation see Ref. [5]. It can, however, get trapped in 
local minima of this non-linear optimization problem, or 
become very slow especially for the dilute systems con- 
sidered here. As an example see the badly converged 
density profiles obtained by standard DMRG approaches 
in Fig. 2. 

In our implementation of the MG-DMRG algorithm, 
we start by constructing the target lattice model and a 
hierarchy of models on coarser grids. Starting from the 
coarsest level we optimize the wave function and inter- 
polate it to the next finer level, repeating this procedure 
until we reach the target system. Many generalizations 
are possible, for example iterating the procedure by going 
back to coarser levels, or starting from the finest instead 
of coarsest level. 

The restriction operation maps a system to a coarser 
grid, merging n (typically n — 2) sites into one. The 
model, given by the Hamiltonian H and defined in the 
local basis {a}, is mapped to a restricted model H in a. 
truncated local basis {a} for the n merged sites. The 
truncation, denoted by an isometry T^^^^^^, is straight- 
forward for continuum models and an approach for lat- 
tice models will be discussed below. Any error due to 
the truncation will be corrected when returning to finer 



scales, as long as we stay in the same phase. As illus- 
trated in Fig. 3 the restriction transforms a matrix prod- 
uct state Ai, . . . , An into 
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Prolongation is the inverse of restriction and maps a 
solution from a coarse grid to a finer one. The isometry 
is inverted and T~^ replaces one index by k new indices. 
From this tensor we can recover a (not unique) MPS 
representation by repeatedly applying singular value de- 
compositions and splitting it into matrices Ai, . . . , An 
(see Fig. 4). It has turned out to be useful to perform a 
standard DMRG update on the newly obtained matrices 
immediately after prolongation while keeping the rest of 
the system on the coarse-grained lattice [' ]. 

As a first application we apply MG-DMRG to bosons 
in a one-dimensional continuum system with an exter- 
nal optical lattice potential, V{x) = Vocos^(fcx), with 
k — n/a and a the size of a unit cell. The continuous- 
space Hamiltonian describing spinless bosons interacting 
through a (5-potential in a system of L unit cells and 
length La is 



H = 



dxtp* {x) 



2m dx^ 



V{x) 



■ip{x) 



dx ip' {x)iip' {x)iilj{x)iip{x) , 



(3) 



where a boson is created at position x with the field 
operator ip^{x), satisfying the usual commutation rela- 
tions. We express energies in units of the recoil en- 

f-2 . 2 

ergy Er — %^. The interaction g is conveniently 
parametrized by the dimensionless coupling 7 — mg/h^n, 
where n is the density. 

In deep optical lattices the low-energy sector of the 
model can be mapped to an effective single band Hub- 
bard model with one site per unit cell. We are, however, 
interested also in weak optical lattices and thus discretize 
the continuum model on a grid with N points per unit 
cell and spacing Ax — a/N. To discretize the Hamil- 
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Figure 3. The restriction of an MPS state reducing it from 
8 to 4 sites. The contraction along the indices /?, (t\ and CT2 
creates a new tensor A^^ ^^^ with a new index (t, but keeping 
the old connections to the neighboring sites (ai and 0.2)- 



tonian (3) on this lattice we replace the Laplacian by a 
second order finite difi^erence approximation and replace 
field operators by lattice annihilation and creation oper- 



ators ■^iJ^x = (i + 1/2) Ax) = 



/Ax 



We end up with a 



Hubbard-likc model in a spatially varying potential: 



H{Ax) = - t{Ax) Y^ [c|c»+i + h.c.j + Y^ ^J.^{Ax)ni 

i i 

^n,(n,-l) (4) 



C/(Ax) 
2 



with t{Ax) = (fi2/2m)/Aa;^ [/(Ax) = .g/Ax, and 
/ij(Ax) = t/((i + l/2)Ax) + 2(fi2/2m)/Ax2. A similar 
lattice model can be formulated for fermions. 

With this definition of the Hamiltonian for arbitrary 
Ax the implementation of MG-DRMG is straightfor- 
ward. The matrix elements of the isometry T are 
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where Ci, a[ and a are particle-number eigenstates and 
we truncate the maximum occupation of a site at A^maxj 
i.e. ai,a G {0, . . . , A^max}- Due to particle-number con- 
servation, T is a iVniax X -^max block-diagonal matrix. 
Note that we start from a coarse-grained lattice and per- 
form only prolongations. 

As a benchmark for the multigrid algorithm wc con- 
sider an optical lattice with Vo/Er — 6 and I/7 = 0.1 
corresponding to the insulating phase [ ] at unit filling. 
Our MG-DMRG simulation were performed with up to 
A^ = 128 lattice sites per unit cell (Ax = 0.0078125) 
with A'lnax = 2 keeping M — 200 states and using 12 
sweeps of single-site updates at each level. We also per- 
form standard DMRG simulations starting from either a 
random initial state, a state obtained from an infinite- 
size growing procedure, or a few steps of imaginary time 
evolution (not shown). The infinite-size growing proce- 
dure is commonly used to obtain good initial states for 
one-dimensional systems and has been proven to be very 
efficient in most cases. We use the implementation of 
the ALPS DMRG code [20], which performs the growing 
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Figure 4. Prolongation operation doubling the number of 
sites: first we transform the basis a of the initial matrix A 
into two new local bases cti and (T2. With a singular value 
decomposition we then split the rank-4 tensor A into two 
matrices Ai and ^2- 
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Figure 5. Comparison of the energies in a bosonic optical 
lattice [Vo/Er = 6, I/7 = 0.1) with L = 32 unit cells dis- 
cretized with increasing discretization A*' = 16, 32, 64, 128 ob- 
tained with different strategies: DMRG with initial state 1 
optimizing an initial random state, DMRG with initial state 
2 initializing the system with an infinite size procedure and 
linearly increasing the number of states (results obtained with 
the ALPS DMRG code [20]), MG-DMRG, and MG-DMRG 
combined with local optimization in the prolongations. 



procedure on a state with very small bond dimension and 
increases the bond dimension linearly with the number 
of sweeps thereafter. 

Fig. 2 shows the density profile obtained with the 
three approaches. Clearly, the standard DMRG ap- 
proaches are trapped in configurations with a glob- 
ally non-homogeneous density distribution and further 
sweeps are not effective in redistributing the particles. 
The multigrid method, on the other hand, achieves a 
symmetric distribution, since it performs the optimiza- 
tion of the large-scale features on the initial coarse mesh 
with just N — 2 sites per unit cell, where convergence 
is very fast. Subsequent calculations on finer lattices are 
initialized with the prolongated solution of the coarser 
lattice. This is already close to the ground state, and only 
the local fine-scale features of the wave function need to 
be optimized. 

The better convergence of the MG-DMRG is also re- 
flected in the energies shown in Fig. 5. While for modest 
discretizations, standard approaches yield energies com- 
parable to MG-DMRG, they encounter severe conver- 
gence problems for smaller values of Ax where the differ- 
ence between multiple scales of the dilute system become 
more and more important. The most reliable method is 
MG-DMRG combined with optimization in the prolon- 
gation. 

MG-DMRG opens new interesting applications for 
DMRG that have not been accessible before. As an ex- 
ample we combine MG-DMRG with time evolution [7- 
], to study heating caused by non-adiabaticity when 
ramping up the amplitude of an optical lattice. We 



start from the ground state of a homogeneous system 
of length L = 16 and A^ = 16 grid points per unit cell 
calculated by MG-DMRG. We evolve it in time using a 
fourth-order Trotter decomposition with At = Q.Olh/Er- 
Non-adiabaticities due to ramping at a finite speed cause 
heating and we plot the energy difference to the ground 
state in Fig. 6 for three different ramp profiles and several 
total ramping times. For the calculation of the ground 
state energies in weak optical lattices MG-DMRG was 
used. We observe that as the ramp speed decreases, dif- 
ferences in ramp shape are less important than the to- 
tal ramping time, indicating that the exact shape of the 
ramp profiles play a minor role in experiments and exper- 
imentalists should focus on determining optimal ramping 
times. 

In DMRG simulations of weakly doped t-J or Hubbard 
ladder models [23-2")] the hole density shows similar con- 
vergence problems as seen above for dilute particle sys- 
tems. In particular it has been observed that for six holes 
in more than 2 x 64 sites the standard DMRG algorithm 
fails to distribute the three bound hole pairs evenly over 
the ladder [26], and MG-DMRG can be of use here. The 
restriction step of mapping the model to a coarser lat- 
tice is, however, not as straightforward as in continuum 
models. We propose to use the contractor renormaliza- 
tion (CORE) method [27] to find a good approximation 
of the model in the reduced Hilbert space of the coarser 
models, and to iterate this procedure in further restric- 
tion steps. For the specific case of doped ladder mod- 
els, the first step maps 2-site rungs or 4-site plaquettes 
to a hardcore boson model for the hole pairs, or an ex- 
tended plaquette model containing hole pairs, magnons. 
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Figure 6. Heating due to finite ramping speeds in a system of 
length L = 16 for I/7 = 0.08 with M = 400. Different Unes 
show the energy difference to the ground state while ramping 
up to Vo/Er = — >■ 4. (inset) Ramping functions used in the 
time evolution: linear, Vo{t)/Er = it/ty (dashed); exponen- 
tial, Vo{t)/Er = 4[exp(4f/iv) - l]/[exp(4) - 1] (dotted); and 
s-like, Voit)/Er = 4 [3{t/tvf - 2{t/tvf] (solid). 
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Figure 7. Coarse-graining a ladder of spin-1/2 fermions into 
a system of hard-core bosons. Spin singlets are mapped to 
empty sites, and hole-pairs to hard-core bosons in the first 
step. Further restriction steps of the bosonic model are similar 
to dilute bosonic models discussed above. 



and holes [2n]. Further restriction steps map to simpler 
bosonic models for the hole pairs, as illustrated in Fig. 
7. After prolongation back to the full lattice model the 
ground state wave function can be further improved by 
repeating the multigrid scheme can be performed. Now 
one can use knowledge of the approximate ground state 
to perform the restrictions of the basis, instead of using 
CORE. Details of this method and results of this ap- 
proach will be published elsewhere. 

We point out that MG-DMRG is a fundamentally dif- 
ferent approach from the one taken by tree-tensor net- 
works [29] or the multi-scale entanglement renormaliza- 
tion ansatz (MERA) [30, 31]. In those approaches, a 
new class of wave functions is proposed that describes 
the system at several levels of coarse-graining and all 
levels are optimized simultaneously - which can still suf- 
fer from convergence problems at fine scales. Instead, 
our approach relies on standard matrix-product states 
which can be optimized and evaluated much more easily 
and much faster, but uses a hierarchical coarse graining 
to achieve a faster and more reliable optimization than 
standard DMRG. 

Our algorithm can be easily combined with other opti- 
mization schemes for DMRG, such as using iTEBD [8, 32] 
to directly simulate the thermodynamic limit. One 
can also easily generalize the restriction (2) and pro- 
longation to tensors of higher rank, in order to apply 
the multigrid scheme to other tensor network states, 
e.g. MERA [30, 31], projected entangled pair states 
(PEPS) [33] and infinite PEPS (iPEPS) [34]. 
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knowledges the ANR Grant No. 09-BLAN-0097-01/2. 
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